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Reduced frequencies. 


^ kM 
U ’ B 


Modified potential difference across the lifting surface 
Reference length 

Direction cosines of a normal to the lifting surface 


Local lift 
Slope of a line 
Mach number 

Unit normal or summation variable 
MspU^ dynamic pressure 

^ij 

Generalized aerodynamic coefficients = 

Surface of integration 

Hyperbolic radius squared (X^-X)^ - (Y^-Y)^ - (Z^-Z)^ 
Dimensional time 

Ut 

Nondimensional time, — ^ 

Air stream velocity 


Induced velocity components, also termed as backwash, sidewash and 
normalwash components 
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x,y,z 

X,Y,Z 
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Influence coefficient matrix relating normal velocity and velocity doublets 
Dimensional space coordinates 

X y z 

Nondimensional space coordinates, — , — , — 

2p i e 

A transformation matrix, Eq. 48. 

Constant of a line 

= yM^- 1 

Displacement normal to the lifting surface = gjj" exp ( i k T) 

= Y„-Y 
= X^-X 
Air density 

Velocity potential, a scalar quantity 
_ 0gik Mx Modified potential 

Circular frequency of harmonic motion, radians per second 


INTRODUCTION 

Accurate determination of unsteady aerodynamic forces is essential for precise evaluation of 
the aeroelastic stability characteristics of flight vehicles. In the subsonic case, the integral formula- 
tion is simple and computational methods have been well developed (ref. 2). In the supersonic case, 
the problem is complicated by the fact that pressure discontinuities arise because of the conical flow 
field emanating from geometric irregularities such as cranked leading edges, wing tips and interfering/ 
interacting surfaces. A closed form solution for such problems is unlikely, and numerical superposi- 
tion methods of finite element type, using sources, velocity potential doublets or pressure doublets 
as basic variables have received attention in the literature (refs. 3-24). 

The source superposition method gives a very simple integral relationship between the po- 
tential and the downwash field (which is determined by the mode shapes) in the non-interacting 
case (ref. 3). However, in the interacting case, the potential is first related to the source strength 
and this is in turn related to the downwash distribution (ref. 4). Thus, by this method, two sets of 
equations are involved in solving the problem. In addition, integration over wake regions and non- 
unique “diaphragms” is necessary as a part of the solution. 

In the velocity potential method, there is a direct relationship between the downwash (the 
mode shapes) and the velocity potential (ref. 1 ). Diaphragm regions are no longer necessary and the 
wake regions do not need detailed modeling since their behavior is determined by the trailing edge 
potentials of the wake-producing surface. The integral relations which are more complicated than 
in the source superposition method have been considerably simplified in the current work. 

The pressure potential or kernel function method is a ‘direct’ approach via a relation between 
downwash and pressures and aerodynamic coefficients (refs. 5-8). This integral relation is, however, 
even more complicated than in the velocity potential approach. 



For arbitrary configurations, the numerical methods employed may broadly be classified as 
collocation and finite element methods. Collocation methods assume, a priori, certain mode shapes 
or series expansions of the unknown parameters such as pressure and doublet strengths. The coeffi- 
cients of these series or modes are determined from a set of algebraic equations established by satis- 
fying the integral relation only at an appropriate number of collocation points. The number of 
equations is comparatively few and is computationally efficient for simple configurations (ref. 9). 

In a recent paper, Cunningham (ref. 10) uses the three-dimensional kernel function with a judicious 
selection of pressure functions for interacting surface configurations. However, the application of 
pressure collocation methods to handle general configurations is complicated by the difficulty of 
choosing pressure modes for complex multi-dimensional shapes such as wing-body combinations. 

In finite element methods, the integration over the dependence domain is replaced by a sum 
of integrations over a number of simple elemental domains (finite elements). Over each area element, 
the unknown parameter is expressed as a sum of simple functions. A number of finite element 
shapes have been used, such as squares, Mach or characteristic boxes, and triangular or quadrilateral 
elements. Numerical approaches differ also in the choice of functional variation within each element 
and in the integration method over the element (refs. 1 1-20). 

In Mach or characteristic box schemes, planform edges have usually been approximated by 
jagged representations which result in erratic behavior of the pressure over the whole surface. More 
recent versions of Mach box programs are described in references 1 1 and 12. 

Stark (ref. 13) used elementary characteristic boxes in developing a digital computer pro- 
gram in which special consideration was given to the handling of subsonic singular leading edges. 
These modifications, however, detract very significantly from the basic simplicity of the Mach or 
characteristic box approach, the computational price paid for additional accuracy being large. 

A triangular representation of the dependence domain using a linear distribution of sources 
was developed in references 14 and 15. This method offers acceptable accuracies with far fewer 
elements than other methods. 

Allen and Sadler (ref. 17), using characteristic elements, developed a doublet superposition 
method based on Jones’ integrated potential formulation for planar configurations. They expressed 
the kernel (sine and cosine) functions as parabolic interpolation functions within each element. 
Woodcock and York (ref. 18) extended this approach to interacting wing and wing-body configura- 
tions. 

The integrated potential approach was further developed in reference 23, using linearly 
varying potential doublets within triangular elements. Closed form integrals were employed to 
evaluate the singular functions, while numerical integration methods were adopted for more com- 
plex but analytic functions. 

Although good results were obtained with fewer elements, an arbitrary wing with a control 
surface could not satisfactorily be idealized. This disadvantage of the velocity potential method led 
to the development of the potential gradient method described in reference 25 and this report. In 
this scheme, the potential gradient in the stream direction is considered as an independent variable 
and is assumed to be constant over an element. This results in fewer integrals than in the velocity 



potential method discussed in reference 23. Once again closed form integrals have been obtained for 
singular functions and recurrence formulae derived for non-singular terms. Two sides of a typical 
quadrilateral element are taken parallel to the stream, and computational efficiency is increased be- 
cause the integrals along these two lines and along the boundary of the area of the wing cut by the 
Mach cone vanish. This type of element has also facilitated the provision of automatic grid genera- 
tion. Velocity potential distributions and generalized aerodynamic coefficients have been obtained 
by the use of the potential gradient method and compared with available results derived by other 
methods. 

GENERAL AERODYNAMIC ANALYSIS 

In the present analysis the coordinates x, y, z and time t are replaced in nondimensional form 
by X, Y, Z and T, respectively, where 



( 1 ) 


1/2 

8 being the standard length, U the airspeed, and /3 = [M^ - 1 ] where M is the Mach number. The 
effect of the above transformation is to change the planform of the wing in such a way that its chord 
is lengthened while its lateral dimensions remain the same. At the same time the Mach lines in the X, 
Y, Z coordinate system are inclined at ±45° to the X axis as indicated in Figure 1. 
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S DENOTES AREA CUT OFF BY 
MACH CONE WITH VERTEX AT P 


Figure 1. Domain of Influence in Supersonic Flow 

Next let us suppose that a delta wing with subsonic leading edges is oscillating in the airstream 

with frequency co rads/sec! Then if the velocity potential of the disturbed flow, is replaced 

by a nondimensional modified potential 4 j such that 


4 


( 2 ) 


* = gikM^X/^ 
UC 


it can be proved (ref. 1 ), that 


02 $ 02 $ 
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02$ 

az^ 


+ k'2 $ = 0, 


(3) 


where k' = and k = Furthermore, it is shown in reference 
(3) may be expressed as the integral relation 


that the solution of Equation 


<i>(Xo,Yo,Z^) = - //k(X,Y) 

where the integral is taken over the part of the wing cut off by the Mach cone with vertex at X^, Y^, 
Zq. The symbol K = -^^)/2ir denotes the difference between the modified velocity potential 

above and below the wing and 


cos k'R 
R 


dXdY 


(4) 


R = [(Xq - X)^ - (Yq - Y)2 - (Zq - Z) 2 ] 


(5) 


The surface of the Mach cone with vertex at Xq, Yq, Zq is defined by R = 0 and the mean position of 
the wing is assumed to be in the plane Z = 0. 

In general, the modes of motion of the wing are assumed to be known. The displacement 
normal to the wing’s surface will be denoted by 17^ [ = Ct? exp (ikT) ] where t? is a function of x and 
y. The condition for tangential flow over the wing may then be expressed as 


^ = 11 eikT 

dt 9n 


or 

dr] 

dT 

where 
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dT 
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and — is the velocity normal to the surface induced by the doublet distribution over the wing and 
on 

the wake. The factor exp (ikT) is cancelled through Equation (2) and subsequently. 



Wherf the wing lies approximately in the plane Z = 0 and the displacement tjj ) is small, the 
above relation simplifies. In terms of the modified potential <E> and the nondimensional coordinates 
defined by Equation (1), the boundary condition to be satisfied at a typical point Xq, Yq, Zq is then 
obtained by differentiating Equation (4). For the case considered, Equation (6) is replaced by 
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axo J 


gik'MX 


o 






cos k'R 
R 


dXdY 


(7) 


When the appropriate K distribution that satisfies the above equation has been determined, the values 
34> 

of and can be deduced from Equation (4) by differentiation. The actual velocity com- 
0c4 3(4 

ponents may then be derived by differentiating Equation (2). 

3Xq 3Yq 

The Modified Upwash 


In order to determine the modified velocity components, it is convenient to express cos k'R 
in series form 


cos k'R 
R 
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Z "^2n 


R 


2n-l 


( 8 ) 


n = o 


where 


(_l)n |,-2n 


*-2n - IHT 

The modified upwash W is then given by 


(9) 


W = 


3«1> 

azi 


Z C2„W 


2n 


( 10 ) 


where 


W2n = -J -2 //K(X,Y)R2n-ldXdY 


3Z 

w 

and n = 0, 1,2, etc. 

When n = 0 and Z = 0, 


(11) 


rr ^ . 

*0 = - ii-2 JIt 


( 12 ) 
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where ^ = X^-X and 77 = Yq-Y. By integrating by parts, it may be deduced that 


W 

“ 3V 


ffW 


where 


L = ylog, 


Let us next assume that the area of integration S is divided into a number of small quadri- 
lateral elements E with chordwise sides parallel to the J or X axis. It is further supposed that 
0K 0K 

( = - — ) is constant over each element. Then, since L = 0 over the part of the boundary of S on 
05 oX 

the surface of the cone R = 0, it can be deduced that is given approximately by 


V' 9K 9 r r S 

^ ^ ^ If 

^ s 


the sum of the contributions from all the E elements. The above equation, after integration with 
respect to ^ then yields 

9K r 3 / ^o^ \ 

*0 • I W * zj) ” 


where ^ denotes the contour integral around a quadrilateral element in the anti-clockwise direction 
as indicated in Figure 2. In Equation (16), ^ is replaced by mr? + a where m and a have the values 
corresponding to the particular side of the quadrilateral over which the integration is being performed 
and 


R = |^(m^ -1)57* + 2m a 77 + J 


It should be noted that m^ -1 can be negative and that ■ When the above expression for R is 

substituted in Equation (16), it can be shown that 


Wo = - E .0 


where I_ is given by 


*0 " [ -(m'-DFo + m^Zo" F, -maF^j 
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Figure 2. Contour Integration Along Discrete Elements 


where 


Fo 

F. =i 


dri 

~R 


dr] 

(t?2+Zq^)R 

T?cl77 

(r?^+Zo^)R 


and ]^. indicates ^ around each element E, Figure 2. 
Since 


^ _ mZg" -T?g 
dr? ~ (tj^+Zo^)R 

it follows that the value of F^ can be derived from the relation 


aF2 = mZg^ F| - L 


( 20 ) 


(21) 


( 22 ) 
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Hence, only F_ and Fj need be evaluated. It can readily be deduced that 




(m^ - 1)^^ R + (m^ - 1) T?+mo: 
- 1 


, . . . m > 1 


1 . _i (1 -m^) Tj -ma 

sin 


. . . m<l 


. . . . m = 1 


The integral Fi can also be derived by substituting 7} = tan (0+7) in the integrand and putting tan 
7 = mZja. After some reduction, the required formula can be shown to be 

1 r • -1 f ' *"^0 

L I (7}^+Zg^) [q:^+Zq^ (m^ -1] 


(^+R) 

(5-R) 


where the integral is taken along the line J = mrj + a from tJi to t?j . It should be noted that the 
values of m and a are different for the lower side of the quadrilateral element. 

When the values of all the F integrals over the upper and lower sides of each element have 
been determined, the total integral 1^ for any element can be evaluated. From Equation (18), the 
upwash due to all the elements may then be obtained by summation. 

Similarly, it can be deduced that 


Wj = — // 

^ 9Z R 

t ax 


where 


V oK 
E 3X 


T7L(m??+2Q:) Rrj Fq 

h(n) = ^ - Y + — [a^ -3(m^ -DZ^^] 


Z 2 F, 

+ (3mZ^ -4a^) - 


ImcxZQ Fj 
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is taken around the contour C of each element of area (Figure 2). The corresponding formula for 
W 2 „, for n > 2 is given by 


V 3K 

W2„= ( 2 n-l)L — l2n 

E 

where 

l2„ (H2„-2-<2n-3)Z„>H2„.4)d, 


in which 

«2n 

with 


j^2n+ 1 
2n(2n+l) 


(2n-l) 

2n 




> ^2n-2 


and 


Ho = - R 



(V+Zq^)Ho 

2 


(27) 

(28) 

(29) 

(30) 

(31) 


The integrand in 12^ is analytic througliout the dependence domain and closed form integration is 
possible. However, it is more economical to evaluate the integral numerically, say, by the method of 
Gaussian quadrature. It is clear from the expressions for H 2 n that the integrals 12^ vanish on the 
Mach boundary R = 0. This eliminates the need Tor the determination of the intersection of the 
dependent domain by the Mach hyperbola. However, for hyperbolic radius R» 1 and n large, the 
magnitudes H 2 n grow in geometric progression. The convergence of W for large values of reduced 
frequency k is then slow. However, for practical values of k and R convergence is obtained with a 
reasonable number of terms. 

The Modified Backwash and Sidewash Components 


may 


The corresponding velocity components ^2n’ ^2n ^^o^g the OX and OY axes, respectively, 
also be deduced. 

It can readily be proved that 


3<ho 


E 
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3K 

3X 

(mFj + «F, ) 

a<i>2 _ 
3X„ 
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II 
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O 
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E 

3K 

lx 

(t?L - mZ^ Fj - aZ^ Fi + aF^) 


~dX^ 


2n 


= U 2 n = (2n-l)Z, 


E 


dr? 


n^2 


(32) 

(33) 


(34) 
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where P 2 „ is defined by 

P2n = f 

^R^n-l ■ (2n-l)(r?^+Z„^)P2n-2 

2n 

witJi 

Po = L 


(35) 

(36) 


and 


fR-(7?2 +Zq^)L 

P 2 = (37) 

Similarly, the sidewash components are given by 


3cI>o 

W X/ — 1 V 

y 

3K 

/ R \ 

(38) 

3Yq ° ° 
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E 

3X 



34>, 
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3K 


(39) 
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E 

3X 

[(mj? + «) L - R] 



and 


3<D 

3Y 


2n 
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V2n = + (2n-l)Zo 
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3K 

3X 


H2n-2-n>2 


(40) 


The Wake Field 


In the wake, since it can sustain no lift. 


3X ' /3 


This implies that 


ijc (X-Xje) 
K(X) = K(Xjg)e ^ 


(41) 


(42) 


where K(Xj£) is the value of K at the trailing edge of the wing section being considered. If the wake 
is assumed to lie in the plane Z = 0, its contribution to the normal modified velocity component can 
readily be deduced since 


ax 


i kK(XTE) 

~1 



(X-Xje) 


(43) 


and KCXyg) is the sum of the values for all the elements upstream. Using the ^ distribution 

given by Equation (43), the contribution to the velocity components at a receiving element from the 
wake elements can be written as. 


W 


2n,w = W2n 


(44) 


k k (45) 

V2„,w = 7 ^“PWj«i-fTE» ''2n 

where = (X^-Xj) is the relative upstream distance between the influencing and 

receiving elements 

and = (Xg-X-pp) is the relative upstream distance between the receiving element 

and center of the trailing edge. 

For n = 0, 1,2, etc. the expressions for W 2 fj and V 2 f, ^ are the same as for the lifting 
surface elements. The summation in Equations (44) and (45) denotes the contribution from all the 
wake elements between the trailing edge and the intersection of the Mach hyperbola with the wake 
sheet (see Figure 3). 

Finally, the modified normal wash distribution using Equation (6) may be expressed in 
matrix form 


,ik'MXo = 
dT 


ax 


+ [W^] 



(46) 


The trailing edge values of the velocity potential Kpp can be approximated by 


Kte = [X] 


ax 


where the transformation matrix X is given by 


(47) 


[X] 


in which 



tr = [AXq, AX, .... AXpJr 


(48) 


(49) 
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Figure 3. Consideration of Wake Sheet 

tj. is a row vector defining the distances between the centers of the elements in the rth span from 
which trailing edge the wake sheet emanates (see Figure 3). 

Using Equation (47), the modified velocity potential gradient can be related to the given 
boundary conditions as; 


[W + . X] 


dK 

ax 


d?7 ik'MXo 
dT ® 


(50) 


Calculation of the Velocity Potential 

3K 

Equation (50) can be solved for the velocity potential gradient It is then desired to 

dX 

determine K and subsequently the velocity potential difference A0 = (j>^ -0g. 

From Equation (2) et seq , 


ax 


1 /dAifi 
27 r\ ax 


+ ik'MAip j e 


ik'MX 


(51) 
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A0 

where A^p = — is the nondimensional velocity potential difference. 
Then 


with 


dAip 

IjT 


+ ik'MAip = \p 


, T 3K , -ik'MX 

* 

Since Aip = 0 at the leading edge, the solution to Equation (52) is given by 


( 52 ) 


(53) 


Aip(X) = e 


■ik'MX Yj 


Chord wise 
Elements 


ik'Mf 

/ «/'(^e ^df 


(54) 


If i// is expressed linearly between two element centers. Equation (54) can be integrated in a closed 
form between element centers, and the contributions summed. Finally, K is determined from 
Equation (2). 


Calculation of the Generalized Forces 


The lift on an element dx dy of the wing’s surface is denoted by £ (x,y) exp (ilcT) and hence 
C (x,y) dx dy = 2^ /32 (X,Y) dX dY ( 55 ) 


where the lift distribution 
£ (X,Y) = 27 t 


aK i k k I 

,ax ' /3 J 


- ik'MX 


(56) 


By the principle of virtual work, the generalized aerodynamic influence coefficient (AlC) 
can be written as 


Qij = 4^ 




0K: i k K 


•J 


Aec 


-ik'MX 


(57) 


where Qjj is the nondimensional form of the coefficients defined by 


Q.. = q£" Q.. 
ij ij 


(58) 


In the above equation, Ag is the area of a quadrilateral element of the wing and the contributions 
from all such elements of both sides of the wing are summed. 
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COMPUTATIONAL METHOD 


This section briefly describes the idealization of the interacting wing surfaces, numerical 
calculation of the velocity components and the solution procedure. More detailed information is 
contained in the user’s and programmer’s manuals (refs. 26 and 27). 

Idealization of Lifting Surfaces 

The choice of the grid system plays an important role in the numerical computation of the 
unsteady aerodynamic coefficients in a supersonic flow field. Triangular, rectangular, quadrilateral 
and Mach characteristic grids have been employed in the supersonic analysis (e.g., references 11, 13, 
14, 17, etc.). While well conditioned and more efficient calculations result from the grids based 
partly on Mach characteristic lines, as shown in Figure 4, the automated generation of such elements 
for a complex configuration such as multiple wings with control surfaces is difficult and proved 
beyond the scope of the current work. 





Figure 4. Typical Characteristic Grid Idealization With Control Surface 


In the present development, trapezoidal elements have been chosen as the basis to model the 
interacting wing configurations with control surfaces. All the elements are assumed to have two 
sides parallel to the free stream. Figure 5 shows the idealization of a typical wing with control sur- 
faces. As previously indicated, this element type has specific advantages for the performance of the 
recurring line imegrals around the perimeter. 



Figure 5. Idealization of Wing Planform for 
Constant Potential Gradient Method 


Numerical Computation of the Integrals 

Integrals of the velocity components given by Equations 9, 25, 27 and 32 througli 40, are 
performed for all influencing elements in the dependence domain, with respect to a receiving point 
at the center of an element. These integrals are taken along the contour of each element. Since the 
sides of the elements are parallel to the streamlines, numerical line integrals are performed only along 
leading and trailing edges of the elements. Furthermore, all the integrals vanish on the Mach cone 
and the need for determining the hyperbolic curves of intersection of the cone with the lifting sur- 
face is avoided. In the case of partial elements, the line integration beyond the Mach hyperbola is 
zero. 

The velocity influence coefficients Wjj 2 ^ which are independent of the reduced frequency k 
are first evaluated for each term in the series expansion. For n = 0 and n = 1 , the integrals are easily 
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1 


expressed in closed form, while for higher order terms, say n >2, closed form expressions exist but are 
increasingly complex. Hence it is more economical to evaluate these integrals by numerical methods. 


At low Mach numbers and for large reduced frequencies, convergence of this series is poor 
for far field elements. Numerical difficulties were also encountered in the evaluation of l 2 ^ integrals 
given by Equation (28) since the magnitudes of the l 2 jj terms became very large. However, for a 
given frequency, the total value of the influence coefficient, i.e.. 


W: 


y 

'J* n = 0, 1,2 . . 

approaches an asymptotic value for far field elements. Figure 6 shows a typical distribution of the 
velocity coefficients for various chordwise and spanwise positions (f, 17) of an influencing element 
i with respect to a receiving element j. 



Figure 6. Asymptotic Nature of Far Field Elements 

For 77 = 0, the asymptotic value is seen to have been obtained for all upstream elements be- 
yond ^>0.4 (region B in Figure 6). A similar asymptotic trend is seen for elements at other span 
stations such as 77 = 0.1 , 0.2, 0.3, etc., with a drastic reduction in magnitude of the velocity com- 
ponents. This behavior is discussed for rectangular elements in the Appendix. 
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A dependence domain C with respect to the receiving point j (Figure 6) can be constrained by 
specifying a certain order of magnitude of the velocity components Wy in relation to Wy such that 
the total solution does not vary by more than, say, 2 percent. Thus, the actual computation can be 
performed only for the region A while the upstream region B is represented by the asymptotic values 
obtained in the region A. In this way increased computational economy has been achieved in the 
development of the computer program (see references 26 and 27). This also eliminates the numerical 
difficulty encountered in the computation of l2„ integrals for the far field elements. 

Solution Procedure 

Accurate determination of unsteady aerodynamic forces on lifting surfaces with controls 
requires that the configuration be defined by a large number of elements which prevents the matrix 
relation given by Equation (50) from being solved economically by inversion techniques. Iterative 
methods have been developed to solve such large order linear systems on digital computers of limited 
memory size. In the present work, an iterative technique developed by Bratkovich and Marshall 
(reference 28) has been chosen as the most efficient one. This method is based on the successive- 
over-relaxation technique, without the need to determine the relaxation factor by trial and error 
runs. 

After the frequency independent coefficients have been determined, the total solution is 
obtained for any frequency within a required range. Computational efficiency can be increased by 
using the previously determined solutions as the starting vectors for different frequencies. 

Reference 27 discusses in greater detail the epu time taken for operations such as: 

( 1 ) determination of the coefficients 

(2) decomposition of the matrix 

(3) solution for each mode, and 

(4) solution for each frequency, etc. 

RESULTS AND DISCUSSIONS 

To assess the solution accuracy and the versatility of the potential gradient method, a number 
of calculations have been performed on an IBM 360/65 computer and compared with available results. 
In all the cases considered here, the wing planforms were represented by trapezoidal finite elements. 
Generalized aerodynamic coefficients and pressure distributions were calculated for various reduced 
frequencies and Mach numbers. The examples include planar and nonplanar configurations with 
control surfaces. 


Isolated Wings 


Double Delta Wing 

The wing planform ( Figure 7a) was represented by 1 50 elements, with 1 5 span stations. The 
pressure distributions along the span for various chord stations are compared with the results of 









ref. 29. The correlation is very good, except at the tip section (see Figure at X = 9.0), where the large 
pressure discontinuity is smeared by the discrete element approach. 

Rectangular and Arrowhead Wings 

The generalized aerodynamic coefficients for rectangular and arrowhead wings (Figure 8a and 
8c) for heave and pitch motions are, in Tables 1 and 2, compared with Fenain’s (ref. 22) and Stark’s 
(ref. 13) methods as reported by Woodcock in ref.21 . The number of elements used in chord and 
span directions are shown in the tables. The generalized aerodynamic coefficients determined for the 
Mach number range M = 1 .04 to M = 2.0 and reduced frequencies k = 0, 0.5, and 1 .0, are seen to be 
in very good agreement, in spite of fewer elements used in the present method. 

AGARD Swept Wing With Control Surface 

Figure 8b shows the planform of the AGARD swept wing with control surface. The wing was 
represented by 147 trapezoidal elements with 15 of them on the Control Surface. The oscillating 
modes of the wing considered in the present calculations were; 

1 . Heave tji =1.0 

2. Pitching about the mid-chord (C/2) t ?2 = X - C/2 

3. Chordwise bending 773 = (X - C/2)^ 

4. Flapping of control surface. 

The generalized aerodynamic coefficients were calculated at M = 1.2, and reduced frequencies 
k = 0.5 and 1.0. 

Table 3 shows the comparison of the present results with those of refs. 13 and 22. In spite of 
the large difference in the number of elements used in the present and the referenced methods, the 
generalized aerodynamic coefficients are in good agreement, except for the loss of accuracy in the 
small order terms. 

When the magnitudes of the real and imaginary parts of a complex number differ significantly, 
it is better to compare two complex numbers using their magnitudes and phases. For example, 
amplitudes of Qj 3 and Q,4 are in good agreement while their phase angles differ only by a fraction 
of a degree. 


Interacting Wings 

Rectangular Wing Folding at 50% Semi Span 

In order to check the accuracy of the out of plane velocity components derived in the 
present approach, a rectangular wing of aspect ratio 4 folded at 50 percent of semi span was con- 
sidered. Lift curve slopes at M =\/2 for various fold angles were calculated, and indicate excellent 
agreement with the theoretical values ref. 30 (Figure 9). Mach box results of ref. 30 shown in this 
figure, overestimate the values of lift curve slopes at low fold angles and underestimate them at high 
fold angles. 



TABLE 1 


AGARD RECTANGULAR WING (PLANFORM FIG. 8 (a)) 
MODES Zj = I.O Z2 = X-C/2 AR = 2.0 


Mach 

No. 

Methods 
(Matrix) or 
(Sp/Ch Pts.)* 

Qii 

0 

II 

k = 

0.3 

k = 

0.6 

Re(Q) lm(Q) 

Re(Q) 

lm(Q) 

Re(Q) 

lm(Q) 


Present (49) 



0.215 

1.070 

0.472 

1.715 


M9 (34/26)** 

1,1 


0.205 

1.060 

0.749 

1.820 


M19 (22/17)*** 



0.189 

1.009 

0.348 

1.639 


Present 


3.978 

3.567 

-0.593 

2.938 

-0.556 


M9 

1,2 

3.951 

3.531 

-0.545 

3.058 

-0.946 

1.2 

M19 


3.750 

3.370 

-0.500 

2.840 

-0.294 


Present 



-0.001 

-0.132 

-0.096 

-0.300 


M9 

2,1 


-0.005 

-0.141 

-0.026 

-0.310 


M19 



-0.005 

-0.131 

-0.107 

-0.285 


Present 


-0.370 

-0.419 

0.157 

-0.403 

0.440 


M9 

2,2 

-0.398 

-0.446 

-0.177 

-0.427 

0.350 


M19 


-0.368 

-0.411 

-0.165 

-0.380 

0.450 


Present (99) 



0.008 

1.128 

0.158 

2.065 


M9 (34/54) •• 



0.034 

1.144 

0.146 

2.083 


M19 (14/22)*** 

1,1 


0.019 

1.088 

0.124 

1.997 


Present 


3.880 

3.955 

0.255 

3.781 

-0.035 


M9 

1,2 

3.787 

4.0 

0.088 

3.806 

0.004 

1.05 

M19 


3.542 

3.80 

0.134 

3.648 

0.036 


Present 



-0.200 

-0.193 

-0.3419 

-0.200 


M9 

2,1 


-0.185 

-0.204 

-0.333 

-0.251 


M19 



-0.174 

-0.197 

-0.320 

-0.244 


Present 


1.339 

-0.544 

0.837 

-0.1616 

0.797 


M9 

2,2 

-1.363 

-0.590 

0.790 

-0.250 

0.825 


M19 


-1.293 

-0.571 

0.749 

-0.247 

0.790 


* Spanwise and chordwise grid points 

** As reported in Ref. 21 using the method of Ref. 22 

*** As reported in Ref. 21 using the method of Ref. 13 



TABLE 2. AGARD ARROW HEAD WING (FIG. 8 (c)) 
MODES: Zj = 1.0, Z2 = X-C/2,AR = 4.0 



Methods 


o 

II 


k = 0.5 

lc = 

1.0 

Mach 

No. 

(Matrix Size 
or Sp/Ch Grid)* 


Re(Q) lm(Q) 

Re(Q) 

lm(Q) 

Re(Q) 

lm(Q) 


Present (43) 



0.046 

0.6514 

0.1.61 

1.239 


M9(34, 15)** 



0.059 

0.6085 

0.198 

1.136 


M19(40, 17)*** 

1, i 


0.055 

0.6122 

0.186 

1.144 

2.0 

Present (43) 

1, 2 

1.334 

1.306 

0.1318 

1.254 

0.292 


M9 (34, 15)** 


1.248 

1.222 

0.0766 

1.164 

0.190 


M19(40, 17)*** 


1.255 

1.228 

0.0851 

1.165 

0.208 


Present (43) 



0.027 

0.236 

0.094 

0.436 


M9 (34, 15) 



0.031 

0.226 

0.105 

0.406 


M19 (40, 17) 

2, 1 


0.029 

0.224 

0.096 

0.404 


Present (43) 


0.489 

0.474 

0.0858 

0.448 

0.188 


M9 (34, 15) 


0.471 

0.457 

0.0595 

0.426 

0.140 


M19 (40, 17) 

2, 2 

0.467 

0.452 

0.0650 

0.419 

0.151 


Present (110) 



0.146 

1.104 




M9 (34, 35) 



0.146 

0.913 




M19 (26, 26) 

1. 1 


0.136 

0.887 



1.25 

Present (110) 


2.182 

2.316 

0.156 




M9 (34, 35) 


2.002 

1.911 

0.138 




M19 (26, 26) 

1, 2 

1.936 

1.852 

0.147 




Present (110) 



0.089 

0.442 




M9 (34, 35) 



0.075 

0.330 




M19 (26, 26) 

2, 1 


0.072 

0.329 




Present (110) 

2, 2 

0.858 

0.949 

0.101 




M9 (34, 35) 


0.758 

0.710 

0.100 




M19 (26, 26) 


0.748 

0.705 

0.102 




Spanwise and chordwise grid points 
** As reported in Ref. 21 using the method of Ref. 22 
*** As reported in Ref. 21 using the method of Ref. 13 


TABLE 3 

AGARD SWEPT BACK WING WITH CONTROL (FIG. 8 (b)) 
MACH NO. = 1.2 TOTAL CPU = 306 SEC (IBM 360^5) 




li 

t-0.5 


k-1.0 


>. j 

Ref 22 

Ref 13 

Present 

Ref 22 

Ref 13 

Present 

No. Elements 

34 x51 

20 X 30 

147 

34 X 51 

20x30 

147 


1, 1 

0.0110 

-0.0228 

0.1758 

-0.3690 

-0.4042 

-0.6150 


2, 1 

-0.2780 

-0.2832 

-0.2332 

-0.7220 

-0.7085 

-0.8567 


3, 1 

-0.1620 

-0.1643 

-0.1193 

-0.3500 

-0.3397 

-0.3840 


4, 1 

-0.0100 

-0.0090 

-0.0094 

-0.0230 

-0.0222 

-0.0186 


1,2 

3.8110 

3.6714 

3.515 

4.1910 

4.0734 

4.0721 


2,2 

0.2870 

0.2442 

0.0906 

0.8580 

0.8016 

0.8924 


3,2 

0.7160 

0.6824 

0.6056 

0.9520 

0.9106 

0.9666 


4, 2 

0.0200 

0.0186 

0.0131 

0.0260 

0.0235 

0.0305 


1,3 

3.5220 

3.3533 

3.916 

2.7140 

2.5672 

2.6444 


2,3 

2.8130 

2.7026 

2.958 

2.2980 

2.2108 

2.0063 

O 

3, 3 

1.444 

1 .3568 

1.547 

1 .2660 

1.1923 

1.1051 

-1 

< 

4,3 

0.0850 

0.0803 

0,0672 

0.0830 

0.0784 

0.0526 

LU 

oc 

1,4 

0.5900 

0.5723 

0.5563 

0.5600 

0.5437 

0.5313 


2,4 

0.5320 

0.5101 

0.4984 

0.5020 

0.4834 

0.4748 


3,4 

0.4830 

0.4583 

0.4494 

0.4550 

0.4332 

0.4244 


4,4 

0.0530 

0.0514 

0.0448 

0.0490 

0.0478 

0.0419 


1, 1 

3.4790 

3.3506 

3.2200 

3.6460 

3.5473 

3.4314 


2, 1 

0.0410 

0.0084 

0.1699 

0.5440 

0.5016 

0.5390 


3, 1 

0.5810 

0.5532 

0.4637 

0.8420 

0.8084 

0.8644 


4, 1 

0.0160 

0.0146 

0.0074 

0.0250 

0.0227 

0.0311 


1,2 

1 .6770 

1 .7325 

1 .2420 

1.6100 

1 .5832 

1 .8967 


2,2 

2.4910 

2.4589 

2.4289 

1 .7880 

1 .7385 

1 .8345 


3, 2 

1.3350 

1 .3055 

1.2182 

0.9040 

0.8632 

0.8562 


4,2 

0.0780 

0.0738 

0.0764 

0.0610 

0.0590 

0.0485 


1,3 

-1 .2090 

-1.1954 

-1 .0223 

0.1750 

0.1563 

-0.0139 


2, 3 

-1.0740 

-1 .0507 

-1.2672 

0.1180 

0.0996 

0.0826 


3,3 

-0.2560 

-0.2489 

-0.3390 

0.4570 

0.4360 

0.5124 

o 

4,3 

0.0220 

0.0198 

-0.0049 

0.0440 

0.0412 

0.0373 


1,4 

-0.0410 

-0.0298 

-0.0170 

-0.0270 

-0.0168 

-0.0054 


2,4 

-0.0320 

-0.0220 

-0.0104 

-0.0180 

-0.0093 

0.0010 


3,4 

-0.0260 

-0.0149 

-0.0173 

-0.0120 

-0.0023 

-0.0059 


4,4 

0.0020 

0.0032 

0.0042 

-0.0040 

0.0052 

0.0058 


Note: Modes are defined in the text. 



(PER RAD) 



Figure 9. Lift Curve Slope For An Aspect Ratio 4.0 Rectangular Wing 
With Folded Tips 


AGARD Wing-Tail Configuration 

In the interacting case, the generalized aerodynamic coefficients were calculated for the 
AGARD wing-tail coplanar configuration (Z = 0) shown in Figure 10. Four antisymmetric modes 


of the form: 

Wing 

Tail 

Mode 

T?1 

Y(X-2.25|Y|-0.85) 

0 

wing twist 

T?2 

Y|Y| 

0 

wing bending 

7?3 


Y 

tail roll 

1?4 


1Y| (X-3.35) 

tail pitch 

were considered. 





The wing and tail were represented by 1 1 1 and 71 trapezoidal elements respectively. 
Generalized aerodynamic coefficients were calculated at M = 3.0 for reduced frequencies k = 0 and 
1.5. Table 4 compares these results with those of ref. 24 and 3 1 . The generalized aerodynamic 
loads on the wing due to wing modes, and on the tail due to tail modes are in excellent agreement 
with the referenced methods. However, the interference effects, i.e., loads on tail due to wing modes, 
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TABLE 4 

AGARD WING-TAIL INTERFERENCE M = 3.0 Z = 0.0 
NO. ELEMENTS =111 (WING) + 71 (TAIL) TOTAL CPU = 246 SEC (IBM 360/65) 


■ 

k = 0.0 

k = 

1.5 


■ 

U 

Ref 24 

Ref 31 

Present 

Ref 24 

Ref 31 

Present 


1., 1 

-0.0226 

0.0189 

-0.0187 

0.0966 

0.1066 

0.0901 


2, 1 

0.3035 

0.2789 

0.3287 

0.3846 

0.3238 

0.3895 


3,1 

-0.2152 

0.2226 

-0.1075 

-0.0394 

0.1438 

-0.0695 


4, 1 

-0.1550 

-0.0006 

-0.0843 

-0.0147 

0.1438 

-0.0360 


1, 2 




-0.0700 

-0.0668 

-0.0819 


2, 2 




-0.0759 

-0.0530 

-0.0888 


3,2 




-0.1531 

0.1701 

-0.0400 


4, 2 




-0.1033 

0.0216 

-0.0179 

d" 

1, 3 
2,3 








3, 3 




0.0168 

0.0127 

0.0140 

Ul 

cc 

4,3 




0.0050 

0.0008 

0.0032 


1,4 








2, 4 








3,4 

0.4650 

0.4338 

0.4756 

0.4517 

0.3859 

0.4541 


4,4 

0.2882 

0.3018 

0.2904 

0.2965 

0.2578 

0.2903 


1, 1 




0.1486 

0.1345 

0.1593 


2, 1 




0.0890 

0.0865 

0.1083 


3, 1 




0.0769 

-0.0612 

0.0007 


4, 1 




0.0559 

-0.0612 

0.0008 


1.2 




0.0309 

0.0463 

0.0300 


2, 2 




0.2363 

0.2040 

0.2498 


3,2 , 




0.0239 

0.0670 

-0.0003 


4, 2 




0.0197 

0.0398 

0.0030 

1,3 







a 

s 

2, 3 
3,3 




0.2560 

0.2283 

0.2635 


4, 3 




0.1786 

0.1669 

0.1820 


1,4 








2,4 








3,4 




0.1632 

0.1518 

0.1820 


4, 4 




0.2188 

0.1910 

0.2300 


Note: Modes are defined in the text. 



do not agree in magnitude among the methods compared. The discrepancy may be resolved by 
conducting additional correlation with other methods. 


Tail-Fin Configuration 

A T-tail configuration, with tail mounted on fin a Z3 = 1 .2 as shown in Figure 10 was considered 
as the next example with the following modes of oscillation: 



Tail 

Fin 

Mode 

ni 

= 0, 

Z^ 

Fin bending 

ri2 

= 0, 

Z(X-0.875Z-3.0) 

Fin twist 

V3 

= Y, 

0 

Tail roll 


Tail and fin were represented by 42 and 63 trapezoidal elements. Generalized aerodynamic coeffi- 
cients were calculated at M = 1 .6 for reduced frequencies k = 0 and 1 .5, and are given in Table 5. It 
was intended to compare with the results of reference 32. For some reason Table Q3.1 1 of reference 
32 was not available in that report. However, the corresponding result with fin reflection is shown 
in Table 5. There seems to be no agreement in any of the generalized aerodynamic coefficients. The 
results of ref. 32 appear to be 2 to 3 times higher than the values obtained from the present method 
with no fin reflection. 


TABLE 5 

AGARD TAIL - FIN INTERFERENCE M = 1 .6 Z = 1 .2 
NO. ELEMENTS = 42 (TAIL) -t- 63 (FIN) TOTAL CPU = 93 SEC (IBM 360/65) 




k = 

0.0 

k = 

1.5 


i.i 

Ref. 31* 

Present 

Ref. 31* 

Present 


1, 1 



1 .3389 

-0.0125 


2, 1 



0.0664 

-0.0605 


3, 1 



0.4630 

0.056 


1, 2 

2.7068 

0.8089 

1.4164 

0.6516 

a 

2, 2 

0.3511 

0.0970 

0.2823 

0.1031 

"m 

3, 2 

0.5621 

0.2258 

0.1392 

-0.1359 


1, 3 



0.1251 

0.0370 


2, 3 



0.0211 

0.1795 


3, 3 



0.0350 

0.0025 


1, 1 



4,7124 

0.5200 


2, 1 



0.1324 

0.0449 


3, 1 



1.6933 

0.1209 







d" 

1, 2 



-0.4249 

0.1430 

0 

2, 2 



0.1817 

0.1643 

s 

3, 2 



-0.2333 

0.0047 


1, 3 



1,6239 

0.0527 


2, 3 



0.0074 

0.0166 


3, 3 




0.6700 

0.4457 


*The results of Reference 31 were obtained with fin reflection. 
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Wing-Tail-Fin Configuration 

This was the last in the series of the examples considered. The tail was placed at Z = 0. Four 

modes; 

i) wing twist 

ii) tail pitch 

hi) fin bending, and 

iv) fin twist, 

as defined in the previous examples were considered. Generalized aerodynamic coefficients were 

calculated at M = 3.0 for frequencies k = 0 and 1 .5 and are presented in Table 6. For this configuration, 
no other results are available for comparison. 


TABLE 6 

AGARD WING - TAIL - FIN INTERFERENCE M = 3.0 Z = 0 
NO. ELEMENTS = 73 (WING) + 70 (TAIL)+ 1 10 (FIN) TOTAL CPU = 203 SEC (IBM 360/65) 



k = 

0.0 

k = 

1.5 

i, j 

Real (Qjj) 

IMG (Qjj)/k 

Real (Qjj) 

IMG (Qjj)/k 

1, 1 

-0.0202 


0.0839 

0.1641 

2, 1 

-0.1333 


-0.1779 

0.1444 

3, 1 

-0.0980 


0.0492 

0.0066 

4, 1 

-0.0460 


0.0240 

0.0185 

1, 2 





2, 2 

0.3420 


0.3137 

0.2160 

3, 2 

-0.01756 


0.0042 

0.0084 

4, 2 

-0.0414 


-0.0059 

0.0170 

1, 3 





2, 3 



-0.0046 

-0.0006 

3, 3 



-0.0420 

0.2918 

4, 3 



-0.0217 

0.0361 

1, 4 





2, 4 

-0.0144 


-0.001 7 

0.0059 

3,4 

0.3400 


0.03314 

0.0780 

4,4 

0.0545 


0.0551 

0.0711 


Influence of Element Mismatch on Pressure Distribution 

The lifting surfaces are commonly represented by discrete elements bounded by continuous 
lines drawn from root chord to tip chord and leading edge to trailing edge. For highly tapered wings, 
this may result in an undesirable concentration of small elements at the tip. Elements of nearly 
equal areas could be obtained by terminating some lines at partial span stations (e.g., see the control 
surface in Figure 5). 
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In subsonic doublet methods such a discontinuity creates undesirable fluctuations in the 
pressure distribution (ref. 2). Some misgivings were voiced concerning similar effects in the present 
method. They do not seem to be significant as evidenced by the following examples. 

Figure 1 1 a and b shows the chordwise pressure distribution at two spanwise stations adjacent 
to the leading edge crank of the double delta wing of Figure 7. Two cases of ‘element mismatch’ are 
shown. In both the cases the pressure distributions are relatively smooth. Figure 1 lb also shows the 
pressure distribution with element downwash points located at 0.75 and 0.8 of the element chords. 
No significant effect was observed. However, in the case of wings such as Figure 8c, pressure dis- 
tributions with the downwash point taken at 0.5 of the element chord gave unsatisfactory results. 

In general, a choice of downwash point anywhere between 0.6 to 0.8 of the element center chord 
should give satisfactory results. 



CONCLUSIONS 

The potential gradient approach, unlike the pressure potentials, results in simplified integral 
formulae. By expanding the kernel, the singular integrals have been evaluated in closed form without 
the need for principal or finite part integral techniques. Integrals vanish on the Mach cone and the 
need for determining the hyperbolic curves of intersection of the cone with the lifting surface is 
avoided. Furthermore, idealization of lifting surfaces by trapezoidal finite elements with two sides 
parallel to the streamlines, requires numerical integrals to be performed along the other two sides 
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only. Such elements are amenable to automatic grid generation schemes. Comparable results to other 

methods are obtained with far fewer elements. 
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APPENDIX 

BEHAVIOR OF THE FAR FIELD ELEMENTS 


In Figure 6, it appears that the velocity influence coefficient has an asymptotic trend for 
increasing f- 

Suppose in Figure a rectangular sending element has sides (2e, 2y) and coordinates ($, rj) 
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For T}>0, the significant terms are like 

„ k'" 

k CnS+— • • 

16 

and there is no obvious reason why an asymptotic behavior occurs. This should receive more detailed 
study. 
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